/**
  CNOK project, Anyang Normal University, IMP-CAS

  \class TAInterpolate
  \brief Polynomial interpolation. Note that this is a mathematical tool class.
  \author SUN Yazhou, asia.rabbit@163.com
  \since 2020/07/09
  \date Last modified: 2020/09/03 by SUN Yazhou
  \copyright 2020-2023 SUN Yazhou
  \copyright CNOK project, Anyang Normal University, IMP-CAS
*/

#include <cmath>
#include "TAMath.h"
#include "TAException.h"

/// Polynomial interpolation using Neville's algorithm
/// given n points in arrays x and y, and input x, this routine returns the
/// interpolated func value y at xx, and assigning the error estimate to dy.
/// Ref. Numerical Recipes in C: p109
template <typename T1, typename T>
T1 TAInterpolate<T1, T>::PolyInter(const T *x, const T1 *y, int n, const T &xx, T1 *dy){
  // find the element in array x that is closest to xx
  int nm = 0;
  T dx, dxm = fabs(xx-x[0]);
  T1 c[n], d[n]; // c: Pi..(i+m)-Pi(i+m-1), d: Pi..(i+m)-P(i+1)..(i+m)
  for(int i = n; i--;){
    if((dx = fabs(xx-x[i])) < dxm){
      dxm = dx; nm = i;
    } // end for and if
    c[i] = y[i]; d[i] = y[i];
  } // end for over i
  if(fabs(dxm/xx) < 1.e-14) return y[nm]; // xx equals to x[nm], return y[nm]
  ///// then the Neville's algorithm /////
  // nm-- so as to accommodate Pi..(i+m)=P(i+1)..(i+m)+dm,i
  T1 result = y[nm--], ddy; // the initial approximation to y; ddy: the correction
  // for m=0, c[i+1]-d[i] reduces to y[i+1]-y[i]
  // so it is ok to initalize c and d to y. c and d starts from m-1
  for(int m = 1; m < n; m++){ // loop over level m
    for(int i = 0; i < n-m; i++){
      T dxi = x[i] - xx, dxmi = x[i+m] - xx;
      T1 den = dxi - dxmi; // x[i]-x[i+m]
      // two input x's are (to within roundoff) identical
      if(0. == fabs(den)) TAException::Error("TAInterpolate",
        "PolyInter: Elements in array x too close to each other.");
      den = (c[i+1] - d[i]) / den;
      // update c and d to level m
      c[i] = dxi * den;
      d[i] = dxmi * den;
    } // end for over i
    // note that the following line is to choose the path, i.e. whether c or d at
    // this turning point. nm+1 is the position of xx. If it is above the center
    // of the current P column (2*(nm_1)<n-m), we need to downwards, i.e., choosing
    // c, otherwise choosing d to go upwards.
    // As choosing c means going downwards, the relative position of nm does not
    // change, so nm+1 is the original nm as we have made nm-- in the beginning
    // while choosing d means going upwards, the relative position of nm decreases
    // by 1, so use nm is just, and decrement nm as the same time, hence nm--
    result += ddy = 2*(nm+1) < n-m ? c[nm+1] : d[nm--];
  } // end for over m
  if(dy) *dy = ddy; // the error estimator

  return result;
} // end member function PolyInter

// Polynomial interpolation, same as the above, but using Lagrange's algorithm
template <typename T1, typename T>
T1 TAInterpolate<T1, T>::PolyInterLag(const T *x, const T1 *y, int n, const T &xx, T1 *dy){
  T r = 0., d, m; // the final result, the denominator and the numerator
  for(int i = 0; i < n; i++){
    d = m = 1.;
    for(int j = 0; j < n; j++){
      if(i == j) continue;
      if(fabs(x[i]-x[j])/x[i] <= 1.e-14) TAException::Error("TAInterpolate",
        "PolyInterLag: Elements in array x too close to each other.");
      m *= xx   - x[j];
      d *= x[i] - x[j];
    } // end for over j
    r += m/d*y[i];
  } // end for over i
  if(dy) *dy = 0.; // the error estimator
  return r;
} // end member function PolyInterLag

/// \param len is the length of array x (or y), so that the program would choose
/// the closest interval to envelope xx in the center
/// \param np has the same meaning as in the other overload of PolyInter
/// \param warn when used for extrapolation, warn set to false can supress out-of-bounds warning
template <typename T1, typename T>
T1 TAInterpolate<T1, T>::PolyInter(const T *x, const T1 *y, int len, int np, const T &xx, T1 *dy, bool warn){
  if(np > len) TAException::Error("TAInterpolate", "PolyInter: required points > length of the input arrays");
  int i = TAMath::locate(len, x, xx, np, warn); // the proper subsript to start with
  return TAInterpolate<T1, T>::PolyInter(x+i, y+i, np, xx, dy);
  // return TAInterpolate<T1, T>::PolyInterLag(x+i, y+i, np, xx, dy);
} // end member function PolyInter
